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ABSTRACT 

As a result of resonance overlap, planetary systems can exhibit chaotic motion. Planetary chaos 
has been studied extensively in the Hamiltonian framework, however, the presence of chaotic motion 
in systems where dissipative effects are important, has not been thoroughly investigated. Here, we 
study the onset of stochastic motion in presence of dissipation, in the context of classical perturbation 
theory, and show that planetary systems approach chaos via a period-doubling route as dissipation is 
gradually reduced. Furthermore, we demonstrate that chaotic strange attractors can exist in mildly 
damped systems. The results presented here are of interest for understanding the early dynamical 
evolution of chaotic planetary systems, as they may have transitioned to chaos from a quasi-periodic 
state, dominated by dissipative interactions with the birth nebula. 



1. INTRODUCTION 

The presence of chaotic motion in planetary systems 
is well established. As in numerous other dynamical sys- 
tems, chaos in pla netary orbits a ppears as a result of 
resonance overlap ( Chirikov |1979 ) . For small bodies in 
the solar system, clustering of va rious mean-mo tion reso- 
nances leads to chaotic diffusion ( Wisdom|1980 ). Indeed, 
consequences of chaotic motion can be observed both in 
the asteroid belt, as well as Kuiper belt. The motion 
of the planets in the outer solar system is also chaotic, 
but is well bounded, so th e system is stable over ex - 
tremely long periods of time (Murray and Holman|1999|. 



Conversely, in the inner solar system, secular resonances 
drive chaotic motion, and the excursion s in orbital ele- 
ments can be quite large (Laskar 1989). In particular, 



it has been shown that Mercury s proximity to the 
secular resonance may lead to a dramatic i ncrease in its 
eccentricity, followed by eventual ejection fLaskarJll9"96| 
Batygin and Laughlin 2008 Laskar and Gastmeau[2 009 ) . 



All planets in the solar system reside on orbits that 
are relatively far away from the sun, and as a result, 
form a nearly undissipative Hamiltonian system. As 
the discoveries of extra-solar planets have mounted, it 
has become apparent that a large class of planets re- 
side in close proximity to their host stars. Similarly 
to the Galilean sattelites system, tidal dissipation plays 
an important role in the dynamics of these, "hot" ex- 
oplanets. As a direct consequence of tidal dissipation, 
motion of close-in exoplanets has been assumed to be 
regular (Wu and Goldreich 2002). The same can be 
said for planets and dust particles whose orbital eccen- 
tricities and inclinations are constantly damped during 
early epochs of planet formation. In particular, the pres- 
ence of gas gives rise to dissipation i n the form of stokes 



drag ( Beauge and Ferraz-Mello|1993 ). Additionally, non- 
uniform reemission of absorbed sunlight gives rise to dis- 
sipati ve Poynting-Roberts on drag for ~ ^m-sized par- 
ticles (Gonczi et al. 1982). Thus, little effort has been 
directed towards the study of chaos, outside of our solar 
system. In particular, investigations of chaotic motion 
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in planetary orbits, in presence of dissipation, remains 
a sparsely addressed problem. In this study, we seek 
to bridge this gap, with an eye towards identification of 
the dynamical "route" that planetary systems take be- 
tween highly dissipated regular motion and chaotic mo- 
tion within a Hamiltonian framework. It is important 
to note that this examination has direct astrophysical 
implications for understanding the dynamical evolution 
of planetary systems which transition to chaos from a 
quasi-periodic state that is dominated by dissipative in- 
teractions with the nebula, as the gas is slowly removed. 

Classical examples of transition from regular to chaotic 
motion can be found in the context of simple dynami- 
cal systems, such as the Logistic Map, which is usually 
applied to population dynamics, and the Duffing Oscil- 
lator, which describes the motion of a forced pendulum 
in a non-linear potential. In both of the mentioned ex- 
amples, chaos is approached via the "period doubling" 
route, although it is noteworthy that other bifur cations 
that lead to chaoti c motion exist (see for example Albers 
and Sprott (2006) and the references therein). In the 



context of the period doubling approach to chaos, as the 
degree of dissipation is decreased the periodic orbit, char- 
acterized by a period P, onto which the system collapses, 
suddenly changes into a new periodic orbit of period 2P. 
When this happens, the periodic orbit transforms into 
one with two loops, infinitcsimally close to each other 
and to the original shape of the orbit. However, as dis- 
sipation is decreased, the twice-periodic nature of the 
orbit becomes progressively more apparent. If dissipa- 
tion is decreased further, at some point, the system dou- 
bles its period again to 4P, and so on. As this process 
is repeated, the period approaches infinity, which is the 
essence of chaotic motion. In the intermediate regime be- 
tween a chaotic sea and a NP limit cycle (where N is not 
too large), resides a dynamically rich structure, known 
as the "strange attractor", which is a fractal, possibly 
chaotic object of intermediate dimensionality. 

In this work, we show that much like the simple ex- 
amples mentioned above, planetary systems also can ap- 
proach chaos via the period doubling route. Further- 
more, we show that in the context of planetary motion, 
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a strange attractor can exist, given realistic parameter 
choices which loosely resemble that of the solar system. 
Our approach to the problem lies in the spirit of classical 
perturbation theory, where orbit-averaging is employed 
and only a few relevant terms are retained in the Hamil- 
tonian. In addition to yielding deeper insight into the 
physical processes at play, this approach is necessary for 
an efficient exploration of parameter space. Indeed, de- 
spite the considerable advances in computational tech- 
nology in the recent decades, direct numerical integra- 
tion of dissipative systems remains considerably slower 
than that of Hamiltonian systems (because symplectic 
mappings cannot be used), rendering parameter explo- 
ration a computationally expensive venture. We begin 
with a brief review of integrable, linear secular theory for 
a coplanar planetary system, accounting for eccentricity 
dissipation in section 2. In section 3, we extend our anal- 
ysis to non-linear secular perturbations and demonstrate 
the appearance of global chaos in a purely Hamiltonian 
framework. In section 4, we add dissipative effects and 
show the period doubling approach to chaos, and the ex- 
istence of the strange attractor. We discuss our results 
and conclude in section 5. 



2. LINEAR SECULAR THEORY 

Consider the orbit-averaged motion of a test particle, 
forced by an eccentric, precessing, exterior planet. We 
can envision the orbital precession of the planet to be 
a consequence of perturbations from yet another com- 
panion^), which is too distant to have an appreciable 
effect on the test particle under consideration. For the 
sake of simplicity, let us fix the precession rate, g, and 
the eccentricity, e p , of the perturbing planet to be con- 
stant in time. If the test particle is far away from any 
mean-motion commensurability with its perturber, we 
can write its secular Hamiltonian as 



e sin Azu 



X -j]{h 2 + k 2 ) + ^/3(h 4 + k A ) + 1 {hh p + kk p ) 

w 

where n is the test particle's mean motion, a is its semi- 
major axis, h — esintu and k — e cos w are the compo- 
nents of the eccentricity vect ors, and the subscript p de - 
notes the perturbing planet ( Murray and Holman||1999| . 
In the above Hamiltonian, r), p and 7 are coefficients that 
depend on masses and semi-major axes only, and their 
functional forms are presented in the appendix. In the 
regime where 77 does not overwhelmingly exceed other pa- 
rameters, a Hamiltonian of this form is often re ferred to 



as the second funda mental model for resonance ( Henrard 
and Lemaitre|198 3), and also describes mean-motion res 



onances in the planetary context, although the variables 
take on different meanings. As will be discussed below, 
the fourth-order term introduces non-linearity into the 
equations of motion and renders the m non-integrable 



allow ing for the appearance of chaos (Lithwick and Wu 
20T0I). 



"Although the purpose of this study is to investigate 
the onset of chaos, it is useful to first consider the reg- 
ular, integrable approximation. Thus, let us neglect the 
non-linear term (i.e. set /3 = 0) for the moment. An ap- 
plication of the linear form of the perturbation equations 




Fig. 1. — Phase space portrait of a test particle, forced by an 
exterior m = 15m^ perturber, in the linear, integrable approx- 
imation. The blue curves depict un-dissipated orbits, while the 
opaque gray line shows a dissipated orbit, with 8 = 0.02?). The red 
dot onto which the dissipated orbit converges represents the fixed 
point, which acts as a global attractor for the dissipated system. 



to the Hamiltonian, yields the equations of motion. 



dh _ 1 dH 
dt na 2 dk 



dk _ I dH 
dt na 2 dh 



(2) 



Let us now add dissipation into the problem. In plan- 
etary systems, dissipation may come about in a number 
of ways, but is most commonly discussed in the con- 
text of tidal friction and interactions of newly formed 
bodies with a gaseous nebula. Both of these processes 
lead to a decay of eccentricity and semi-major axes. In 
the case of tides, the semi-major axes decay time-scale 
usually greatly exc eeds that of the eccentrici ty, since 
r a = a/a — e 2 T circ ( Murray and Dermott||l999 ). Conse- 
quently, the decay 01 semi-major axes can be neglected in 
most circumstances. The same is generally true for the 



dissipative effects of the nebula ( Lee and Pcalc 2002 1 , al 
though the formalism may be somewhat more complex. 
Consequently, we model the damping of the eccentricity 
as an exponential decay with a constant circularization 
timescale: de/dt = —Se, where 5 = l/r circ , while we 
neglect the decay of semi-major axes altogether. Intro- 
ducing complex Poincare variables, z — eexpi-cu, where 
i = v 7 — T) equations (2) with the inclusion of t he dissi- 
pative term can be written in a compact form (Wu and 
Goldreich|2002|>: 



dz 



irjz + i^e p e 



igt 



(3) 



This equation of motion admits a stationary periodic so- 
lution z = jz p / (g — r\ — i8), which can be expressed as a 
fixed point in terms of the variable z = z/z p . Note that 
in z, the system (1) becomes autonomous. Physically, 
this fixed point corresponds to a state where the eccen- 
tricity of the particle is constant, while its apsidal line is 
co-linear and co-precessing with the perturbing planet. 
Whether the particle is apsidally aligned or anti-aligned 
with the planet depends on the sign of (g — rj) . 

The secular fixed point has been discussed in some 
detail in exoplanet literature. It has been shown that 
co-planar systems can approach the fixed point on 
timescales considerably smaller than that of typical plan- 
etary system lifetimes, given sufficient tidal dissipation 
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Fig. 2. — Phase space portraits of a test particle, forced by an exterior m = 15m^ perturber, as given by non-linear secular theory. The 
corresponding perturber eccentricity for each portrait is labeled. Red dots indicate fixed points of a given portrait. The thick green curve, 
where present, depicts the separatrix. The opaque gray arrows mark the basins of attraction of each stable fixed point. 



| |Wu and Goldreich|[2002| |Mardling||2007| ). The resulting 
disappearance of one of the secular eigen-modes from the 
dynamics of a multi-planet system can yiel d dynamical 



stabil ity, where it is otherwise unachievable ( Lovis et al. 
20111. The significantly non-zero inner eccentricity of 



a fixed point has been invoke d to explain ongoing tidal 
dis sipati on in close-in planets ( |Mardling|2007| |Batygin"et 
|al.||2009[ ). Furthermore, for close-in planets, where tidal 
precession and general relativity play dominant roles, the 
exact value of the fixed point eccentricity proves to be 
a function of the planetary Love number and planetary 
mass. This allows one to infer infor mation about a tran- 
siting planet's interior from its orbit ( |Batygin et al.|2009| 
and resolve the sin(z) degeneracy in non-transiting sys- 



tems (Batygin and Laughlin 20111. 

Here, we neglect general relativistic, rotational and 
tidal precession. This yields perturbation equations that 
are scale-free (i.e. only dependent on the semi-major 
axis ratios), but the particular examples shown here are 
not directly applicable to close-in planets (although the 
extension of the framework to account for additional pre- 
cession is very simple - see ( |Batygin and Laughlin|[2011[ ) 
and the references therein). 

In the parameter regime described, the general solution 
to the equation of motion is 



e P 7 



e (ig+5—iri)t 



(4) 



4 



Batygin & Morbidelli 



where c is an integration constant that depends on the 
initial conditions. In absence of dissipation, the phase 
space portrait is a familiar set of concentric curves that 
close onto themselves. However, if dissipation is intro- 
duced in the system, the phase space area occupied by 
the orbit begins to contract. Given a sufficient amount 
of time, the particle settles onto the co-precessing fixed 
point. This is an important distinction between Hamilto- 
nian and dissipa tive systems: Ha miltonian flows cannot 
have attractors ( Morbidelli 2002 1 . The existence of at- 
tractors requires the presence of dissipation. 

Figure 1 illustrates a phase space portrait of un- 
dissipated, as well as damped motion of a test-particle, 
perturbed by an exterior, m = I5m^ planet, orbiting a 
Sun-like (M* = 1M ) star. Variables are plotted such 
that the radial distance depicts the eccentricity of the 
test-particle, while the polar angle represents the angle 
between the apsidal lines of the particle and the planet. 
The blue curves depict un-dissipated orbits, while the 
opaque gray line shows a dissipated orbit, with 8 = 0.02^. 
The red dot onto which the dissipated orbit converges 
represents the fixed point, which acts as a global attrac- 
tor for the dissipated system. The semi-major axes ratio 
between the test-particle and the planet is chosen to be 
a = a/a,p = 1/2, e p = 0.1 and g = 23"/year. 

Recall that the position of the fixed point is a func- 
tion of the perturbing planet's eccentricity. As will be 
apparent below, this is central to our argument. If the 
perturbing planet resided on a circular orbit, the fixed 
point would be at the origin. Furthermore, in our for- 
mulation, whether the fixed point is apsidally aligned (to 
the right of the origin) or anti-aligned (to the left of the 
origin) depends on the precession rate assigned to the 
perturbing planet. 

3. NONLINEAR SECULAR THEORY AND THE ONSET OF 
CONSERVATIVE CHAOS 

Now consider the evolution of the test-particle without 
omitting the fourth-order terms in equation (1). The 
equation of motion now reads 



dz 
dt 



a/1 - \z 2 \ (irjz + i(i\z 2 \z + i 7 e p e ,9i ) - Sz (5) 



Note that with /3 = 0, and the square root expanded 
to first order in e, we recover equation (3). Here, the 
square root appears because we no longer limit ourselves 
to the linear form of Lagrange's planetary equations. The 
purpose of the square root is to correct for the fact that 
(h, k) variables are only a low-eccentricity approximation 
to the true canonical variables, although its inclusion is 
not instrumental to our results. No general analytical 
solution for this equation exists, and one must resort 
to numerical integration to explore the dynamics. As 
before, it is useful to begin the analysis in absence of 
dissipative effects. 

The addition of the non-linear term introduces impor- 
tant qualitative differences into the solution. First and 
foremost, if the perturbing planet is eccentric, there are 
now up to three no n-trivial fixed points presen t, instead 
of one (e.g. Ch.8 of Murray and Dermott| ( |1999| ). One of 
these fixed points can be unstable (saddle point) and re- 
sides on a critical curve (i.e. separatrix) that surrounds, 
both a librating as well as circulating orbits. Figure 2 
shows the phase space portraits of the particle motion, 



perturbed by a planet of the same parameters as before, 
but with different eccentricities. 

If the perturber's orbit is circular (Figure 2A), the sit- 
uation is quite similar to the linear case. In fact, if we 
omit the square root in equation (5), then a simple ana- 
lytical solution exists. In this case, the fixed point is at 
the origin. In direct analogy with the results of the pre- 
vious section, in presence of dissipation, the fixed point 
would attract all orbits. There also exists an eccentricity 
value for the test particle which sets its precession equal 
to that of the perturber. This set of stationary points is 
illustrated in Figure 2A as a red circle. However, as long 
as the perturber's orbit is circular, these stationary con- 
figurations are qualitatively no different than any other 
eccentric orbit. If we now make the perturber slightly ec- 
centric (e p — 0.005), the dynamics changes dramatically 
(Figure 2B). The first new feature is that the circle of sta- 
tionary points breaks in two individual fixed points: one 
unstable at Axu = it and one stable at Atu = 0. A criti- 
cal curve (i.e. the separatrix, green bold curve in panels 
B-E) is generated at the unstable equilibrium point and 
encircles the stable one. Second, the stable fixed point 
that was at the center of the figure moves slightly to the 
left. 

The appearance of new fixed points has ramifications 
for dissipative dynamics. As in the linear example, if 
dissipation (assumed to be finite but much too small to 
noticeably modify the dynamical portrait, i.e. lim <5 — > 0) 
were to be introduced, both of the stable fixed points 
would act as attractors, with their respective basins of 
attraction (shown as gray arrows in Figure 2) separated 
by the critical curve. The stability of fixed points that 
do not lie on the critical curve, can be understood in the 
following qualitative manner. Consider a small libration 
cycle, centered on one of the fixed points. The cycle's 
intersections with the x-axis are placed symmetrically, 
relative to the fixed point. The role of dissipation at the 
higher eccentricity intersection is to decrease the radius 
of libration, while that at the lower eccentricity intersec- 
tion is to increase the radius of libration. Of the two 
antagonist effects, the first wins, because e oc e. Thus, 
the fixed points centered on libration cycles are stable 
foci. 

As the eccentricity of the perturber is increased fur- 
ther to e p — 0.05 (Figure 2C) and then to e p = 0.1 
(Figure 2D) the phase-space area engulfed by the inner 
branch of the separatrix (i.e. orbits centered around the 
stable anti-aligned fixed point) shrinks. Simultaneously, 
the phase-space area occupied by orbits that are librat- 
ing around the aligned fixed point grows. When the per- 
turber eccentricity reaches e p = 0.12, the apsidally anti- 
aligned fixed points collapse onto a single, unstable fixed 
point (Figure 2E) . This implies that if dissipation was to 
be increased, no apsidally anti-aligned attractor would 
exist. If the eccentricity of the perturber is enhanced be- 
yond e p > 0.12, the anti-aligned fixed point disappears 
completely from the portrait (Figure 2D). 

In both "end-member" scenarios we considered (e p = 
and e p — 0.2), in presence of dissipation, only a single at- 
tractor would exist. However, the attractors in these two 
cases arise from different fixed points (i.e. one can not 
be transformed into another by a change in e p ), centered 
around different branches of the separatrix. Recall that 
the sole fixed point that was present in the e p = case, 
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Fig. 3. — Laplace-Lagrange secular solution of Jupiter, Saturn, 
and a m = 15m,0 perturber at a = 1.85 AU. The black and red 
lines, which never go above e = 0.1 show the eccentricities of 
Jupiter and Saturn respectively. The blue and black curves, which 
attain high eccentricity, represent the exact Laplace-Lagrange so- 
lution and the approximate solution, given by equation (6) for the 
perturber 's eccentricity. 

disappeared, when e p — 0.12. Similarly, The fixed point 
that is present in the e p = 0.2 portrait is not present 
when the perturber's orbit is circular. This has impor- 
tant implications for the motion of the particle when ec- 
centricity of the perturber is not maintained at a constant 
value. 

Consider a scenario where no dissipation is applied, 
but the eccentricity of the perturber is varied adiabat- 
ically between e p — and e p = 0.2. Here, "adia- 
batically" means that the oscillation period of the per- 
turber's eccentricity greatly exceeds the apsidal circula- 
tion/libration period of the test-particle. Regardless of 
the particle's starting condition, its orbit will eventually 
encounter the separatrix. Since the separatrix is an or- 
bit with an infinite preriod, its crossing necessa rily leads 



to chaotic motion (Bruhwiler and Cary 1989). In fact 
the situation is analogous to the motion or an amplitude- 
modulated pendulum. It has been shown that the chaotic 
region of such a system occupie s the phase-space area 
that is swept by the separatrix (|Henrard and Henrarcl] 
1991||Henrard and Morbidelli||1993p . As a result, by en- 
suring that the eccentricity of the perturber reaches zero 
and extends above e p = 0.12 at every oscillation, we en- 
force the entire phase-space within e < 0.6 to be swept 
by the critical curve, causing all test-particles within this 
e-limit to become chaotic. 

Large variations in the perturber's orbital eccentric- 
ity can be induced by secular interactions with a distant 
pair of planets. Consider placing the perturber described 
above at a = 1.85AU, initially on a circular orbit, in pres- 
ence of Jupiter and Saturn, whose initial condition s cor- 
respond to their a c tual o rbits in 1983 (see Ch. 7 of |Mur-| 
ray and Dermott (1999)). The orbital evolution ot the 
massive system can be computed using Laplace-Lagrange 
secular theory, and the resulting solution is presented in 
Figure 3. Note that this solution is approximate at large 
eccentricity, since high order terms are neglected. Due to 
the perturber's proximity to the secular resonance, its 
orbital eccentricity undergoes excursions between e p = 
and e p ss 0.2 on a r sa 2n/(u p — u s ) ~ 1 Myr time-scale, 
where it's are the corresponding eigen-frequencies of the 
perturbation matrix and s refers to Saturn. Furthermore, 
the perturber's longitude of perihelion in this solution 




Fig. 4. — Poincare surface of section, illustrating the chaotic dy- 
namics of the test particle. A section is taken at minimum per- 
turber eccentricity. The blue points correspond to an evolution, 
where the eccentricity of the perturber is given by equation (6). 
Orange points correspond to an evolution, where the eccentricity 
of the perturber varies in a similar manner to that described by 
equation (6), but between e p = 0.05 and e p = 0.2. In the scenario 
where the perturber's eccentricity does not reach zero, the sepa- 
ratrix fails to sweep the entire phase-space, so a resonant island, 
roughly outlined by a black orbit, appears. 

precesses at a nearly constant rate of g = 23" /year. For 
our purposes we approximate the variation in the per- 
turber's eccentricity as 



e p s» 0.2 1 sin(- 



L )*l 



(6) 



Addition of Jupiter and Saturn to the system does not 
significantly modify the evolution of the test-particle be- 
cause of the substantial orbital separation between them 



(«3 



0.2 and a s ~ 0.1). In fact, evaluation of the 



corresponding constants 77, /? and 7 shows that the par- 
ticle's interactions with Saturn can be neglected all to- 
gether, as they only contribute at the ~ 1% level, while 
for Jupiter, it suffices to account only for the additional 
apsidal precession, to which it contributes at the ~ 30% 
level, compared to the effect of the considered planet 
(15m®). Quantitatively, this corresponds to an enhance- 
ment of the coefficients ry and /3, but not 7. As stated 
in the Appendix, where the expressions for the constants 
are given, we have been implicitly retaining the apsi- 
dal contribution due to Jupiter since the beginning of 
the paper, for consistency of the phase-space portraits. 
The difference in longitude of perihelia between the test- 
particle and Jupiter forms a comparatively fast angle, 
and thus can be averaged out. Consequently, we avoid 
its introduction into the Hamiltonian. This is further 
warranted, as the magnitude of the interaction term be- 
tween Jupiter and the test-particle (i.e. 7-,) is about an 
order of magnitude smaller than that of the test-particle 
and the m = 15m® perturbeiR 

Since the introduction of the variation of the per- 
turber's eccentricity, we are now faced with a one-and- 
a-half degrees of freedom Hamiltonian. The dynamics of 
such a system is best visualized by using a Poincare sur- 

1 We could have generated an identical chaotic region by intro- 
ducing another perturber that precesses slightly slower, and tuning 



its parameters, such that the interaction coeff ic ients in the Hamil- 
tonian , 
12010) 

pendulum rather than an amplitude-modulated one. 
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m, 7 are equal (see SidlichovskyJ Il990[l , |Lithwick and Wu| 
3)). Such a system would constitute a frequency-modulated 
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Fig. 5. — This figure depicts the evolution of the fixed point 
and the subsequent approach to chaos. The black curve shows the 
movement of the fixed point on a Poincare surface of section. As S is 
reduced, the fixed point leaves the origin and travels outwards in a 
spiral manner. In the region —5.26 > log ir) 6 > —5.37, a temporary 
2P limit cycle is encountered. The period doubling cascade and 
the onset of chaos (boxed) begins for dissipation rates lower than 
logio 5 = ~ 6 - 4 - 

face of section. As in the phase-space portraits above, on 
a Poincare surface of section, a periodic orbit will appear 
as a point, or a finite sequence of points. A quasi-periodic 
orbit will appear as a curve that closes upon itself, while 
a chaotic orbit will appear as a sea of points, which fill 
a portion of the phase-space. Here, we take a section 
through phase-space every time e p goes through zero i.e. 
with a period of approximately lMyr. 

The Poincare surface of section illustrating the chaotic 
dynamics of the particle is shown in Figure 4. The blue 
points correspond to an integration of the system de- 
scribed above, where the eccentricity of the perturbing 
planet varies according to equation (6). From Figure 4, 
it is immediately apparent that in this setup, the particle 
stochastically explores a large fraction of the phase space 
and no holes appear to exist. We can confirm the chaotic 
nature of this system by measuring its Lyapunov coeffi- 
cient, A, which is a measure of the exponential divergence 
rate of nearby orbits: 



A 



lim 

jV->-o< 



N 

1=1 



\n(di/d ) 
NAt 



(7) 



where do is the initial phase-space separation, di is the 
phase-space separation after some time At and N is 
the number of renormalizations, where the separation 
between t he orbits is manually returned to the initial 
value, do ( |Benettin et al. 1976). Adopting At to be 
the time between successive sections, N = 500 and 
do = 10 -6 , we obtained a positive Lyapunov coefficient 
of A = 4.28 x 10~ 6 years -1 signifying chaotic motion, 
with an e- folding timescale of r ~ (g — if) /A ~ 7 secular 
cycles. Variation of parameters in equation (7) did not 
change our estimates significantly. 

For illustrative purposes, we also performed an inte- 
gration where the eccentricity of the perturber varies in 
a similar manner to that described by equation (6), but 
between e p = 0.05 and e p = 0.2. As already discussed 
above, in such a scenario, we expect that the particle will 
not explore the entire phase-space, as the separatrix will 
fail to sweep all space. The results of this integration are 
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Fig. 6. — Poincare surface of section, showing the period doubling 
cascade and the approach to chaos. As dissipation is progressively 
reduced, the fixed point (organge dot), splits into a IP limit cy- 
cle (gray dots), which subsequently splits into a 4P limit cycle 
(red dots) and finally into an 8P limit cycle (blue dots). Chaos is 
achieved shortly after. 

plotted on Figure 4 as orange dots. In accord with the 
expectations, in this setup, the separatrix falls short of 
sweeping a considerable section of phase-space and as a 
result, there exist islands of stability, which the particle 
never visits. The primary island of stability is outlined 
by a black orbit in Figure 4 and is centered around the 
apsidal libration fixed point of the e p = 0.05 phase-space 
portrait (see Figure 2B). 

4. ROUTE TO CHAOS IN PRESENCE OF DISSIPATION 

Having constructed a system which exhibits chaotic 
motion in the previous section, we can now begin to ex- 
plore the effects of dissipation on chaotic motion. In- 
tuitively, we can expect that in the regime where dis- 
sipation dominates all other effects, no chaotic motion 
can exist. However, the behavior of the orbits in the 
regime that is intermediate between global chaos and 
dissipation-dominated motion, is not apparent a-priori. 

Our strategy is to begin in the dissipation-dominated 
regime and track the behavior of the system, while re- 
ducing 5 in equation (5). The eccentricity evolution of 
the perturbing planet is taken to be governed by equa- 
tion (6). We begin with log 10 <5 = —4. Numerically, this 
regime is one where the orbital cirularization timescale, 
t c = <5 _1 exceeds the free precession rate, 77, by a small 
amount. The solution in this case always falls to a fixed 
point at the origin. This starting point is optimal, since 
increasing dissipation further does not change the loca- 
tion of the fixed point onto which the solution collapses. 
As dissipation is decreased, the location of the fixed point 
begins to depart from the origin in a spiral manner. This 
is shown in Figure (5). 

When dissipation is decreased to log 10 S — —5.26, the 
period of the fixed solution doubles. In other words, the 
solution falls not onto a fixed point, but onto a limit 
cycle. A 2P limit cycle appears as two points rather 
than one on a Poincare surface of section. The series 
of limit cycles, shown as blue points, corresponding to 
—5.26 > log 10 6 > —5.37 is labeled accordingly on Fig- 
ure 5. When the dissipation is lowered below log 10 5 > 
—5.37, the orbit once again collapses onto a fixed point. 
Such behavior is common among systems that approach 
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chaos via period doubling. The trajectory continues to 
collapse onto a fixed point until log 10 8 = —6.4, when the 
period doubles again, and further decrease in the magni- 
tude of dissipation leads to chaotic motion. 

Period doubling is shown in grater detail in Figure (6), 
which is a zoom-in of the box in Figure 5, labeled "ap- 
proach to chaos." At log 10 S = —6.3, the orbit still re- 
sides on a fixed point, shown as a large orange dot. At 
log 10 6 = —6.5, the period has doubled and the limit cy- 
cle is shown as two gray dots. Decreasing the dissipation 
further to log 10 5 = —6.6, each of the two gray dots splits 
into two points, giving rise to a 4P limit cycle. This limit 
cycle is illustrated as four red dots in Figure (6). When 
dissipation is decreased to log 10 S = —6.617, each of the 
four red points splits further into two, resulting in an 
8P limit cycle. This is shown as series of small blue 
points in Figure (6). Decreasing the dissipation process 
repeats the period doubling process. It is noteworthy 
that once the period doubling process begins, the period 
of the limit cycle onto which the solution collapses is a 
very steep function of 5. In other words, the period ap- 
proaches infinity quickly below log 10 S = —6.6. 

As already mentioned in the beginning of the paper, 
an important feature that dissipative systems can ex- 
hibit, which Hamiltonian systems cannot, is the strange 
attractor. In our setup, the strange attractor appears at 
log 10 S — —6.7. The phase-space portrait of a strange 
attractor with log 10 S = —6.8 is illustrated in Figure (7). 
The strange attractor is in a sense an intermediate state 
between a limit cycle and global chaos. Although the 
motion on the attractor itself is chaotic, as can be read- 
ily inferred from comparing Figures (4) and (7), it does 
not occupy the entire available phase space area. This 
is because the attractor has a diminished dimensionality. 
Let us consider the dimensionality of the attractors we 
have encountered thus far. 

Globally chaotic motion, shown in Figure (4) fills the 
entire available phase space area, and thus its surface of 
section lies on a two-dimensional manifold. When strong 
dissipation was introduced into the problem, the motion 
collapsed onto a fixed point which is zero-dimensional 
object. Limit-cycles have surfaces of section of dimen- 
sionality between and 1. A dimensionality of unity is 
achieved if the motion collapses onto a limit-torus, whose 
surface of section appears as a curve that closes upon 
itself. Further decrease in dissipation results in the ap- 
pearance of strange attractors, whose surfaces of section 
lie on manifolds of intermediate dimensionality, between 
1 and 2. 

The dimensionality of the phase-space portrait can be 
related directly to the Lyapunov exponent and thus pres- 
ence of chaos. If motion that originates from different 
initial conditions converges onto a single fixed point or 
limit-cycle attractor, the Lyapunov exponent must be 
negative, signaling periodic motion. For example, in our 
system, a fixed point with log 10 5 = —6 is characterized 
by A = —8.1 x 10~ 7 years -1 . The sign of the Lyapunov 
exponent changes if the dimensionality of the attractor 
exceeds unity. Indeed, motion on the strange attractor, 
shown in Figure (7) is characterized by A = 1.87 x 10~ 6 
years -1 . Note that although motion is chaotic, the e- 
folding timescale corresponds to ~ 16 secular cycles, a 
factor of ~ 2 longer than that of the undissipated sys- 



e sin Aro 




Fig. 7. — A strange attractor. The shown object corresponds to a 
dissipation level of log 10 8 = —6.8 and is characterized by a positive 
Lyapunov coefficient of A = f .87 X 10 — 6 , signaling chaotic motion 
on the attractor. The Minkowski-Bouligand dimensionality of this 
particular attractor is D = 1 .75±0.005. Its persistence requires the 
lack of islands of stability in the occupied region of phase-space. 



tern. 

There are many ways to define a fractal dimension. 
Here, we shall work in ter ms of the Mi nkowski-Bouligand 
dimension (see for example |Ott| ( |l993| ). The Minkowski- 
Bouligand dimensionality of a particular strange attrac- 
tor can be computed by utilizing a box- counting algo- 
rithm. In this approach, the phase-space is divided into 
an even number of sub-regions (boxes) and the number of 
boxes, occupied by the attractor is counted. This is per- 
formed over a large number of scales, typically decreasing 
the box size by a factor of 2 upon each iteration. The 
slope of the line which describes the number of occupied 
boxes as a function of the box size in log — log space is 
the dimensionality of the object. We have performed this 
calculation for the strange attractor presented in Figure 
(7), covering 6 scales. As a result, we find that the at- 
tractor has a dimensionality of D = 1.75 ± 0.005. 

As can be expected from simpler examples, such as the 
Duffing Oscillator, the strange attractor is only present 
for a limited range of parameters. We have performed 
additional experiments where dissipation was decreased 
further. We found that the majority of the attractor 
breaks up into large chaotic regions by log 10 S = —7.7, 
and at log 10 5 = —8, the surface of section is once again 
essentially filled, signaling a return to the global chaotic 
sea, characteristic of the undissipated system (Figure 4). 

In the analysis above, we started from a configuration, 
where in absence of dissipation, chaos engulfed all avail- 
able phase-space, in which no holes existed. Consider 
what would happen if we were to redo the experiment 
using a slightly different setup, that yields a small island 
of stability around the apsidally aligned fixed point, as 
shown in Figure 4. In such a system, the strange attrac- 
tor would not persist, for eventually, the particle would 
necessarily end up in the neighborhood of the island of 
stability, and driven by dissipative effects, find itself in 
the basin of attraction of the fixed point. This implies 
that in configurations where chaos is not global, the pres- 
ence of dissipation tends to guide the constituents to- 
wards regular orbits, with the degree of dissipation di- 
rectly dictating the time-scale needed for the system's 
arrival to a quasi-periodic state. 
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5. DISCUSSION 

In this paper, we investigate the onset of chaotic mo- 
tion in planetary systems, where dissipative effects play 
an important role. Using a semi-analytical perturbation 
approach, we have shown that planetary chaos appears 
through a period doubling cascade. We have further 
demonstrated that strange attractors can exist in the 
context of a planetary problem under the condition of 
global chaos in absence of dissipation. 

It is important to consider the astrophysical signifi- 
cance of this process, beyond purely academic interest. 
As already mentioned in the beginning of the manuscript, 
one can expect the period doubling cascade to occur dur- 
ing early epochs of a planetary system's dynamical evo- 
lution, as the evaporation of the birth nebula leads to a 
gradual decrease in dissipation. As a result, the work pre- 
sented here describes how the dynamical portrait of a sys- 
tem may evolve shortly after formation, when gas is grad- 
ually taken away. Note that the dissipation timescale, 
corresponding to the strange attractor, is typical of a 
late-stage proto-plan etary disk (i.e. r circ ~ 10 6 years) 
( Lee and Peale|[2002 ) . In other words, the example con- 
figuration described here may correspond to a planetary 
system where the planets are already massive enough to 
have essentially decoupled from the gas, and are forcing 
a small planetesimal, which still feels considerable drag. 

Other forms of dissipation, such as tides, are abun- 
dantly present in the planetary context, and are of special 
importance for hot exoplanets. This is further relevant, 
since understanding the dynamics of multiple close-in 
planets is becoming increasingly important, as their num- 
bers in the observed aggregate grow. Particular interest 
is exhibited towards orbital configurations that converge 
to a fixed point, since a stationary system is required for 
obtaining an esti mate for the Love nu mber, /c2, of extra- 
solar gas giants (Batygin et al. 20091. Although linear 
theory predicts that a dissipated system's arrival to a sta- 
tion ary configuration is only a m atter or time (i.e. tidal 
Q) ( |Goldreich and Soter 19661, non-linear theory pre- 
sented here, suggests that one should exercise caution, 
as a fixed point is not always the end-state. 

The same problem can also be turned around. We have 
shown here that limit cycles reside in limited parame- 
ter regimes. Thus, an observed system, whose orbital 



evolution follows a limit-cycle, can be used to place des- 
perately needed constraints on the tidal quality factor, 
which remains among the most unconstrained parame- 
ters in planetary science and whose physical origin is an 
area of ongoing research ( |Wu||2005| ) . 

For close-in planets, the effect of general relativity and 
tidal precession plays in favor of approach to the fixed 
point, rather than any other attractors. This is because 
it enhances the coefficients r\ and j3, but not 7, in the 
Hamiltonian. Since the orbital precession of a putative 
external perturber will generally be comparatively slow, 
the enhanced precession of close-in planets will tend to 
de-tune any resonance. As the relative amplitude of the 
external perturbations is diminished, the dynamics ap- 
proaches the e p = phase-space portrait seen in Figure 
2A. Naturally, this will also lead to stabilization of the 
orbits. It is noteworthy that such an effect is also at 
play in the solar system, as additional precession from 
GR places the free precession of Mercury further from 
the z/ 5 secular resonance, diminishing its chances of e.jec 
tion (Batygin and Laughlin 2008 Laskar and Gastineau 



2009) 



Finally, it is worthwhile to consider the limitations 
of the presented model. Indeed, we have approached 
the problem by utilizing a classical perturbation theory, 
where only a few relevant terms are retained in the dis- 
turbing function. As already mentioned above, this ap- 
proach was necessary for the exploration of parameter 
space, as the efficiency offered by conventional direct in- 
tegration is not sufficient. Although the approach we 
take here breaks down at high eccentricities, in the so- 
lar system, it has been successful in capturing the im- 
portant physical pr ocesses that govern c haotic motion 
(iSidlicho vsky 1990). The recent work of Lithwick and 
Wu| ( |2010p has further confirmed this to be true in the 
case ot Mercury's orbit. Thus, we expect that inclusion 
of the full disturbing function will only modify our find- 
ings on a quantitative level. However, future numerical 
confirmation and re-eavluation of the work done here will 
surely be a fruitful venture, especially if performed in the 
context of a particular observed system. 
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APPENDIX 
Coefficients of the Hamiltonian 

In this work, we choose to write the coefficien ts, such that they appear in t he equations of motion without pre-factors. 
The notation used here is identical to that of Murray and Dermott (1999), i.e. b denotes a Laplace coefficient, a is 
the semi-major axis ratio, and T> = d/ da. Throughout the paper, we account for the induced precession that arises 
from the m = 15m^, a p — 1/2 perturber as well as Jupiter, with aj = 0.178, but not Saturn. Evaluation of the 



formulae below shows that Saturn's effect is negligible, 
the m = 15to0 perturber. The resulting formulae read: 



The eccentricity forcing, taken into account is solely due to 



M* a p 



(2a p V + a 2 p V 2 ) b { °\a p ) + 



mj a 
M* a j 



(2ajV + a 2 jV 2 ) bf{a.j) 



(1) 



n 1 m p a 
32 \M^a^ 



" " ^3 + a 4 p4)6 W K) + ^^ (4a 3 p3 + a 4 p 4 )6 (0) (aj) 



M+aj 

n rrir) a , „ o^o\,ti)> 



(2) 
(3) 
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The first and the second terms in r\ and /3 arise from the m = 15m0 perturber and Jupiter respectively. All terms in 
7 correspond to the m = 15to0 perturber. 
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